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PREFACE 


The  Institute  for  Defense  Analyses  was  requested  to  assist  the  Office  of  Naval 
Research  in  the  planning  and  execution  of  the  Infrared  Analysis,  Modeling,  and 
Measurements  Program  (IRAMMP).  This  document  summarizes  a  portion  of  the  work 
performed  with  ONR  under  Task  A- 180,  "Infrared  Clutter  Characterization  and  Modeling," 
on  alternative  algorithms  for  Infrared  Search  and  Track  (IRST)  systems  during  the  period 
June  1994  to  September  1995.  The  work  was  performed  in  coordination  with 
Mr.  Douglas  N.  Crowder,  Naval  Surface  Warfare  Center  (NSWC),  for  the  Sensor 
Technology  Office,  Special  Projects  Office,  ARPA,  under  the  technical  cognizance  of 
Mr.  Thomas  Wiener. 

This  document  has  not  been  subjected  to  formal  IDA  review. 
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1.  INTRODUCTION 


From  the  Earth's  surface  to  altitudes  that  are  not  too  high,  the  temperature  may 
decrease  with  increasing  height.  If  the  temperature  decreases  at  a  large  enough  rate,  the 
index  of  refraction  will  increase  with  increasing  height.^  Such  behavior  can  cause  a 
mirage,  i.e.,  more  than  one  image  of  an  object  located  at  certain  distances  from  the  viewer. 

This  particular  version  of  the  effect  is  called  an  "inferior"  mirage  because  the  normal 
image  (one  due  to  a  bundle  rays,  each  of  which  connects  an  object  point  to  the  viewing  eye 
without  crossing  another  ray  in  the  bundle)  appears  to  be  above  a  second,  abnormal, 
image.  The  phenomenon  occurs  when  some  rays  emanating  from  an  object  toward  the 
Earth’s  surface  either  begin  with  or  achieve  negative  curvature  before  reaching  the  surface, 
so  that  they  bend  upward,  being  totally  reflected  by  the  atmosphere.  These  rays  then  cross 
one  another,  forming  a  caustic,  before  reaching  the  viewer.  The  result  is  an  inverted  image 
of  the  object  appearing  below  the  normal  image  created  by  other  rays  from  the  same  object 
points.  Each  ray  contributing  to  the  inferior  mirage  image  passes  through  a  minimum 
height  between  the  object  and  the  eye,  but  no  rays  contributing  to  the  normal  image  have 
minima  between  the  object  and  the  eye. 

Another  version,  called  a  "superior"  mirage,  exists  when  an  abnormal  image 
appears  above  the  normal  one  and  is  the  result  of  rays  that  are  totally  reflected  from 
atmospheric  layers  at  higher  altitudes  than  the  viewing  eye.  The  superior  mirage  is  due  to  a 
temperature  inversion  that  occurs  when  the  Earth's  surface,  which  in  this  case  is  often 
water  or  ice,  is  colder  than  the  air  above  it,  and  the  atmospheric  temperature  increases  with 
height  instead  of  decreasing.  Atmospheric  refraction  then  causes  light  rays  to  bend  with 
positive  curvature,  i.e.,  in  the  same  sense  as  the  Earth's  surface.  Rays  that  arrive  from  a 
sufficiently  distant  object  converge  on  the  viewing  eye  from  above,  each  attaining  a 
maximum  altitude  between  the  eye  and  the  object. 

If  the  rays  intersect  one  another,  forming  a  caustic,  so  that  they  arrive  at  the  eye  in 
reverse  order  with  respect  to  height,  the  image  that  they  form  of  the  object  will  appear  to  be 


^  However,  the  rate  at  which  the  temperature  decreases,  the  so-called  lapse  rate,  will  normally  grow 
smaller  with  increasing  height,  eventually  becoming  small  enough  for  the  index  of  refraction  to  reach  a 
maximum  and  then  decrease  with  increasing  height  (cf.  Ref.  1,  pp.  77  ff.). 


inverted.  Other  rays,  arriving  at  the  eye  with  smaller  elevation  angles,  will  have  no  maxima 
between  the  eye  and  the  object  and  will  therefore  maintain  their  initial  order  with  respect  to 
height.  The  two  sets  of  rays  will  then  form  an  inverted,  superior  mirage  image  above  a  true 
image  of  the  object. 

As  an  aid  to  understanding  some  of  the  observed  behavior  of  mirages,  this 
document  considers  some  analytical  connections  between  geometrical  properties  of  such 
phenomena  and  physical  properties  of  the  atmosphere.  Section  II  reviews  the  well  known 
relationship  between  the  index  of  refraction  at  optical-infrared  frequencies  and  the  atmo¬ 
spheric  temperature  as  functions  of  altitude.  Section  HI  considers  the  geometrical  proper¬ 
ties  of  rays  in  a  spherically  stratified,  continuous  propagation  medium,  assumed  to  be  a 
good  model  of  the  Earth's  lower  atmosphere.  Section  IV  discusses  consequences  of  this 
analysis  for  the  nature  of  complex  mirages. 


11.  INDEX  OF  REFRACTION  IN  A  SPHERICALLY 
SYMMETRIC  ATMOSPHERE 


It  is  customary  to  deal  with  the  index  of  refraction  n(r)  of  the  atmosphere  in  terms 
of  the  refractive  modulus  (or  N  unit)  defined  by 


N  =  (n-l)xl06 


(1) 


For  optical  wavelengths  greater  than  0.23  |i  and  for  infrared  wavelengths,  N  for  the  atmo¬ 
sphere  is  given  by  Edlen’s  formula  (Ref.  2,  p.  18-7) 


N  =  Nd-Nw  , 


(2) 


where  the  first  term  on  the  right  side  of  (2)  is  for  a  dry  atmosphere  and  is  given  by 


Nh  = 


ao  +• 


ai 


a2 


1- 


f  \2 
f  V  ' 


1- 


P  (To +  15.0) 


(3) 


and  the  second  term,  which  is  the  contribution  due  to  water  vapor,  by 


Nw  = 


Co 


\2 
V  ' 

.Cl 


w 


(4) 
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In  (3)  and  (4)  P  is  the  total  pressure  in  mb,  T  is  the  temperature  in  °K,  Pq  =  1013.25  mb, 
To  =  273.15  °K,  Pw  is  the  partial  pressure  of  water  vapor  in  mb,  and  v  =  lO^A  which  is 
the  wave  number  in  cm“l  for  the  wavelength  X  in  microns.  Also 

ao  =  83.42, 
ai  =  185.08, 
a2  =  4.11, 
bi  =  1.14  X  105, 
b2  =  6.24  X  104, 

CO  =  43.49, 

Cl  =  1.7  X  104. 

Numerical  examples  assuming  a  value  of  25  mb  for  Pw  and  ground  temperatures  on  the 
order  of  300  °K  indicate  that  the  contribution  of  the  water  vapor  term  Nw  to  the  refractive 
modulus  is  well  below  1  percent  of  its  total  value  and  can  therefore  be  neglected. 

It  is  clear  from  (3)  and  (4)  that  the  behavior  of  the  modulus  of  refraction  N  as  a 
function  of  r  is  the  result  of  the  functional  dependence  on  r  of  the  temperature  T,  the  total 
pressure  P,  and  the  water  vapor  partial  pressure  Pw-  If  the  water  vapor  contribution  is 
neglected  the  functions  P(r)  and  T(r)  will  determine  the  atmospheric  variation  of  the  index 
of  refraction,  by  means  of  which  the  geometrical  properties  of  optical  rays  in  the  lower 
atmosphere  can  be  determined. 

The  hydrostatic  relation  between  the  pressure  P  and  atmospheric  density  p  as  func¬ 
tions  of  the  radial  distance  r  in  an  Earth-centered  polar  coordinate  system  is 


dr 


=  -gP  . 


(5) 


where  g  is  the  acceleration  of  gravity.  The  gas  law  for  an  ideal  gas  is 


P  =  RTp 


(6) 


where  T  is  the  temperature  and  R  is  the  gas  constant.  Eliminating  p  from  (5)  and  (6)  leads 
to  the  differential  equation 


for  which  the  solution  is 


g  ]  ^ 
R,  T 

P  =  Poe 


(8) 


where  re  is  the  radius  of  the  Earth  and  Pq  is  the  atmospheric  pressure  at  the  Earth's  surface. 

In  (8)  g  =  9.8  m/sec2.  Also,  the  gas  constant  R  in  erg/°K-mol  is  8.3145  x  10^, 
and  the  molar  weight  of  air  is  .028971  kgm  (cf.  Ref.  3,  p.l7);  therefore,  since  1  Joule  = 
10^  erg,  in  the  MKS  system 


R  =  286.99389  Joule/°K-kgm  . 


(9) 


Neglecting  the  contribution  Nw  due  to  water  vapor,  it  follows  from  (2),  (3),  and  (8) 
that  the  refractive  modulus  is  given  by 


N  = 


78580 

T(r) 


-.03417  I 


dz 


(10) 


Then,  according  to  (1), 


n(r)  =  l  + 


.07858 

T(r) 


-.03417  I 


dz 

T(z) 


(11) 


III.  RAYS  IN  A  SPHERICALLY  SYMMETRIC  MEDIUM 


For  most  purposes  it  can  be  assumed  that  the  Earth  is  a  sphere  and  that  the  index  of 
refraction  profile  of  the  atmosphere  is  continuously  stratified  in  concentric  spherical  layers 
down  to  the  Earth's  surface.  A  ray  can  be  represented  in  polar  coordinates  relative  to  an 
origin  at  the  center  of  the  Earth  by  an  equation  of  the  form 


r  =  r(e)  . 


(12) 


It  is  assumed  that  the  index  of  refraction  of  the  atmosphere  is  a  function  n(r)  relative 
to  this  polar  coordinate  system.  Then  the  function  r(0)  representing  a  ray  satisfies  (cf. 
Ref.  4,  p.  123  ) 


r(e) 

0  -  00  =  a  j 


dp 


ro  pVpV(p)-a^ 


(13) 


The  one  parameter  family  of  rays  determined,  for  different  values  of  a,  by  (13)  all  pass 
through  the  point  (ro,  0o)-  For  a  given  ray  the  corresponding  parameter  a  can  be 
expressed  in  terms  of  the  angle  \|/  between  the  radius  vector  from  the  origin  of  the 
coordinate  system  to  (ro,  0o)  and  the  tangent  to  the  ray  at  that  point  by^ 


a  =  ±ron(ro)sin\i/  • 


so  that  (13)  becomes 


2  This  is  a  result  of  the  relation  cot\|f  = 


_L*. 

ro  de 


(cf.  Ref.  5,  p.  265). 


0=00 


(14) 


e-00  =±ron(ro)sin\i/  }  — j= 

ro  pVp 


dp 

n^(p)-r§n^(ro) 


If  it  is  contributing  to  a  mirage,  a  ray  must  pass  through  a  maximum  or  minimum 
height  above  the  Earth’s  surface;  i.e.,  at  some  angle  0m  the  radial  coordinate  of  a  point  on 
the  ray  must  satisfy  the  relation 


=  0  . 


This  condition  applied  to  (13)  leads  to 


dr 


=  ± - - '^rJ,n2(rm)-r§n^(ro)sin2v  =  0  , 

0=0„  ron(ro)sin\ir 


which  is  equivalent  to 


rmn(rin)  =  ron(ro)sinv|/  • 


(15) 


Since  (15)  implies  that  rmn(rm)  cannot  be  larger  than  ron(ro),  if  tm  is  a  maximum 
for  the  ray,  n(r)  must  decrease  between  ro  and  rm-  If  I’m  is  a  minimum  for  the  ray,  and 
therefore  no  larger  than  ro,  n(r)  must  be  an  increasing  function  of  r  over  some  neighbor¬ 
hood  of  rm  that  may  or  may  not  include  ro.  In  the  rest  of  this  document  it  will  be  assumed 
that  this  is  the  case  for  r  values  within  some  layer  at  the  Earth's  surface.^ 

The  curvature  of  a  curve  expressed  in  polar  coordinates,  e.g.,  the  ray  given  by 
(13),  is  given  by  (cf.  Ref.  5,  p.  291) 


3  An  example,  discussed  in  Ref.  1,  is  the  case  of  the  unstable  layer  that  usually  forms  on  the  ground  after 
daybreak  and  disappears  at  night. 


K  = 


d^r 


(16) 


Using  the  relations 


—  =  ±^— Jin— 

d0  ron(ro)sin\jr 


^r2n2(r)-r§n2(ro)sin^y 


(17) 


and 


d^r  _  dr  d  dr> 
d0^  d0drvd0; 


(18) 


the  ray  curvatiffe  given  by  (14)  and  (16)  becomes 


K  =  - 


ron(ro)dn 
m^(r)  dr 


sin\|/  . 


(19) 


It  follows  from  (19)  that  the  ray  curvature  is  negative  where  the  slope  of  the  index  of 
refraction  is  positive  and  positive  where  the  slope  is  negative. 

Using  (19)  with  (11),  it  is  possible  to  obtain  the  curvature  of  a  ray  at  any  point  in 
terms  of  the  temperature  profile.  A  straightforward  calculation  provides  the  result: 


K(r)  = 


^  [n(r)-l]ron(ro)sin\|/[^ ^ ^ 

Idr 


m^(r) 


(20) 
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It  is  evident  from  (20)  that  at  a  given  altitude  the  curvature  of  a  ray  is  positive  if  the  slope  of 
the  temperature  profile  is  greater  than  -  .03417,  is  zero  if  the  slope  is  equal  to  -  .03417, 
and  is  negative  if  the  slope  is  less  than  -  .03417. 

Suppose  that  the  eye  level  of  the  viewer  is  at  the  radial  distance  ro  from  the  Earth's 
center  and  the  radial  line  is  the  polar  axis  0  =  0.  Also,  suppose  that  a  ray  begins  at  the 
point  (ro,0)  with  the  direction  of  its  tangent  there  specified  by  an  obtuse  zenith  (measured 
from  the  polar  axis)  angle  X]/-  Then  the  ray  will  be  given  by  (14)  and  (15)  in  the  form 


ro 

®  ~  rnin(rn,)  j 


dp 


pVp^n2(p)-r^n2(rn,) 


,o<e<0 


m 


(21) 


for  r  >  rm,  which  is  the  minimum  radial  coordinate  of  the  ray  as  determined  by  (15),  and  by 


dp 


rm  pVp^"^(P)“rmn^(rm) 


,200,  ^0^0 


m  ’ 


(22) 


for  r  <  rq.  In  (22)  0ni  is  the  angular  coordinate  of  the  ray  point  whose  radial  coordinate  is 
rm.'* 


The  astronomical  horizon,  defined  as  the  direction  orthogonal  to  a  plumb  line,  is 


7t 


given  by  \|/  =  — .  A  declination  angle  ([)  from  the  astronomical  horizon  is  therefore 

2  It 

determined  by  =  (|)  +  — ,  so  that  for  a  small  declination 

2 


sin\|/  =  sin 


xZ 

=  cos<t)~l-—  . 
2 


(23) 


^  Of  course,  if  rm  is  smaller  than  the  Earth's  radius,  the  ray  will  end  at  the  Earth’s  surface  and  will  have 
no  physical  minimum. 


An  optical  horizon  is  defined  by  a  ray  whose  minimum  is  at  the  Earth's  surface  and  is 
therefore  determined  by  an  initial  declination  satisfying  (15)  with  rm  set  equal  to  the  Earth's 
radius  re.  An  approximation  to  the  declination  angle  of  a  horizon  ray  follows  from  (23): 


(|)  ~ 


(24) 


where  the  function  N(r)  is  the  modulus  of  refraction. 


Given  that  the  ray  curvature  K(r)  is  zero  at  rx,  at  any  point  (rm,0m)  for  which  the 
radial  coordinate  satisfies  (15),  whether  rm  is  a  minimum  or  a  maximum  will  depend  on 
whether  rm  ^  rx  or  rm  >  rx.  In  particular,  if  y  is  a  right  angle  the  radial  coordinate  rq  of  the 
ray  associated  with  y  will  be  a  minimum  or  a  maximum  at  the  point  (ro,0),  depending  on 
whether  ro  ^  rx  or  ro  >  rx.^ 


It 

If  ro  >  rx,  when  y  >  —  the  ray  will  start  at  (ro,0)  with  positive  curvature  and  will 
approach  the  Earth's  surface  with  the  index  of  refraction  n(r)  increasing  as  r  decreases. 
When  the  radial  coordinate  of  the  ray  decreases  below  rx  at  which  n(r)  reaches  a  maximum, 
the  index  decreases  with  r  until  the  ray  either  reaches  the  Earth's  surface  or  passes  through 
a  point  (rm,0m)  where  the  radial  coordinate  is  a  minimum.  Up  to  this  point  (21)  determines 
the  ray.  As  the  angle  0  increases  beyond  0m,  (22)  determines  the  behavior  of  the  ray  until 
its  radial  coordinate  reaches  a  maximum. 

7C 

If  \|/  <  — ,  instead  of  (21)  the  ray  equation  becomes 
2 


dp 


p^— — —  (^m) 


,o<e<0 


m 


(25) 


for  ro  r  <  rm,  where  at  the  ray  point  (rm,0m),  the  radial  coordinate  is  a  maximum. 
Beyond  this  point  (22)  determines  the  behavior  of  the  ray  as  0  increases  until  the  ray 
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If  tq  =  Tk,  whether  the  ray  has  a  minimum  or  a  maximum  radial  coordinate  at  (ro,0)  will  depend  on  the 
second,  or  perhaps  even  a  higher  order,  derivative  of  n(r). 


reaches  the  Earth's  surface  or  achieves  a  minimum  radial  coordinate.  If  in  (22)  \\f  equals 
the  angle  Xj/m  that  satisfies  the  equation 


sinxi/m  = 


retire ) 

ron(ro) 


n(re) 

n(ro) 


(26) 


then  according  to  (15)  the  ray  will  become  tangent  to  the  Earth's  surface  and  therefore  will 
determine  an  optical  horizon  for  the  viewer. 

To  examine  the  effect  of  atmospheric  refraction  on  images  of  objects  located  at 
distances  much  smaller  than  the  Earth's  radius,  it  is  useful  to  employ  approximate  versions 
of  equations  (21)  and  (22)  expressed  in  Cartesian  coordinates.  The  origin  of  the  Cartesian 
coordinate  system  is  on  the  Earth's  surface,  so  that  the  vertical  (z)  axis  is  coincident  with 
the  0  =  0  line,  and  the  horizontal  (x)  axis  is  tangent  to  the  Earth.  The  (x,z)  coordinates  then 
may  be  related  to  the  (r,6)  by 


X  =  rO  , 

z  =  r-re  .  (27) 


The  ray  equations  given  by  (21)  or  (25)  in  polar  coordinates  take  the  form 


ZQ  (C+fe 


)i/ft 


dC 


(Q-(zm +re )  V^(zn, ) 


(28) 


where  v(z)  =  n(z  +  rg)  and  z^j  =  r^j  -  r^ .  Dropping  terms  in  (28)  that  are  small  compared 
to  the  Earth's  radius  r^  (28)  becomes 


x  =  ±n(zn,)  j 
ZQ 


dC 

Vn^(C)-n^(zm) 


(29) 
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where  n(z)  is  now  used  in  place  of  v(z)  to  represent  the  index  of  refraction  in  a  plane 
stratified  medium.  In  fact,  (29)  is  the  exact  equation  for  a  ray  in  a  plane  stratified  medium.® 

Starting  with  (29),  a  straightforward  calculation  of  the  ray  curvature  leads  to 


i^(zm)dn(z) 
n^(z)  dz 


instead  of  (19). 

It  is  customary  to  modify  the  index  of  refraction  n(z)  to  accommodate  the  fact  that 
although  rays  in  an  atmosphere  with  a  constant  index  of  refraction  n(r)  relative  to  an  Earth- 
centered  polar  coordinate  system  are  straight  lines,  and  a  horizon  ray  from  the  eye  is 
tangent  to  the  Earth's  surface,  this  cannot  be  the  case  for  a  flat  Earth  relative  to  a  Cartesian 
coordinate  system.  The  solution  is  to  change  n(z)  from  a  constant  to  a  variable  index  of 
refraction  such  that  the  corresponding  horizon  ray  will  curve  enough  to  be  tangent  to  the 
flat  Earth,  i.e.,  by  retaining  a  higher  order  term  in  the  approximation  leading  to  (29). 

A  simple  way  to  accomplish  this  is  to  start  by  getting  rid  of  the  integrals  in  (28)  and 
(29).  Differentiating  both  sides  of  (28)  leads  to 


^  ±  ,  (re  +  ZmMzm) 

V  (’’e  +  v(z)2  -  (re  +  f  (z^ ) 


±  ±  ^  . . . 

+  z  ^^2  _  ^2^^^  j  ^  -  Zm V^(zm  )] 
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Cf.  Ref.  6,  p.  38,  specialized  to  a  plane  y  =  0. 


Assuming  that  v(z)  ~  v(zn,),  it  follows  from  (31)  that 


dz 

dx 


Differentiating  both  sides  of  (29)  and  solving  for  n(z)  leads  to 


n(z)  =  n(zmY  +  |^^j  . 


(32) 


(33) 


Substituting  from  (32)  into  (33)  and  assuming  that  n(zm)  =  v(zm)  provides  the  result 


n(z)~v(zm)- 


1  +  2^  +  ^^^ 


1  +  ^ 


=  v(zni)' 


1  +  ^ 


l  +  ^m 


v(zni) 


1-f-^ 

,  Tey 


l  —  IHL 


so  that 


n(z)~v(zm) 


^  z-z  ^ 
1  +  ^^ 


~v(z„)  +  ^^~v(z)  +  ^^ 


(34) 


‘e  ) 


N.B.,  the  correction  to  be  applied  to  the  index  of  refraction  of  a  plane-stratified 
medium  when  using  a  flat  Earth  approximation  is  ray  dependent:  it  depends  on  the  height 
Zm  at  which  the  ray  has  a  maximum  or  minimum.  The  usual  correction 


n(z)  ~  v(z)-t-- 


(35) 


is  valid  only  for  a  horizon  ray,  which  has  a  minimum  at  the  Earth's  surface  where  Zm  =  0. 


The  temperature  profile  depicted  in  Figure  1  as  a  function  of  the  altitude  z  above  the 
Earth's  surface  will  lead  to  an  index  of  refraction  profile  given  by  (11)  in  the  Earth-centered 
polar  coordinate  system.  When  the  flat  Earth  approximation  is  used,  and  the  resulting 
index  of  refraction  is  modified  in  accordance  with  (35),  the  corresponding  index  of 
refraction  in  a  Cartesian  coordinate  system  with  its  origin  on  the  surface  of  the  Earth  has  the 
profile  shown  in  Figure  2. 


Figure  1.  Temperature  Profile 
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IV.  MIRAGE  CONDITIONS 


According  to  (15),  in  the  Earth-centered  polar  coordinate  system,  for  a  ray  passing 
through  the  eye  at  (ro,0)  to  have  maxima  or  minima  at  several  different  heights  ri,  the  pro¬ 
duct  rin(ri)  must  be  the  same  at  each  height.  With  the  corresponding  Cartesian  coordinate 
system  and  the  flat  Earth  approximation,  for  which  (29)  provides  the  associated  ray 
equation,  the  index  of  refraction  n(zi)  must  have  the  same  value  at  each  height  zi  where  the 
ray  maxima  and  minima  occur. 

This  implies  that  n(z)  must  have  a  minimum  or  a  maximum  between  every  two 
adjacent  heights  z,  and  zi+i  where  a  ray  minimum  or  maximum  occurs.  Moreover,  where 
n(r)  or  n(z)  has  a  minimum  or  maximum  it  follows  from  (19)  or  (30)  that  the  ray  curvature 
is  zero.  In  either  case,  as  the  ray  goes  through  the  height  where  such  a  minimum  or  a 
maximum  occurs,  its  curvature  changes  sign. 

In  addition,  it  is  evident  from  (29)  that  if  n(z)  is  decreasing  with  height  at  the  eye 
height  zo,  the  height  Zm  of  a  ray  maximum  or  minimum  closest  to  the  eye  along  the  ray  path 
must  be  greater  than  zq.  Similarly,  if  n(z)  is  increasing  with  height  at  zq,  Zm  must  be  less 
than  ZQ. 

Figures  3  and  4,  along  with  Figure  2,  illustrate  these  facts.  Figure  3,  which  depicts 
an  inferior  mirage,  shows  a  bundle  of  rays  entering  a  viewing  eye  at  a  height  of  200  feet, 
which  is  above  the  height  where  the  associated  index  of  refraction  profile,  shown  in 
Figure  2,  has  a  maximum.  Figure  4,  which  depicts  a  superior  mirage,  shows  a  bundle  of 
rays  entering  a  viewing  eye  at  a  height  of  80  feet,  which  is  below  the  height  where  the 
profile  has  a  maximum. 


y  ft 
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APPENDIX 


# 


# 


This  Appendix  considers  the  numerical  problem  of  calculating  an  atmospheric 
temperature  profile  for  which  the  associated  index  of  refraction  has  a  maximum,  as  well  as 
rays  with  specified  maximum  or  minimum  heights  in  the  medium.  The  approach  is  aimed 
at  using  Mathcad  Plus  5.0  software  for  Windows  as  the  means  of  implementing  the 
calculations. 


The  formulation  is  based  on  assuming  that  a  cubic  polynomial  with  only  one  real 
zero  defines  the  temperature  profile.  Four  parameters,  one  of  which  can  be  identified  as  a 
scale  factor,  define  the  polynomial,  which  has  the  form 


(A-1) 


where  a,  b,  c  are  arbitrary  parameters  and  sf  is  the  scale  factor.  The  associated  index  of 
refraction  n(z)  is  given  by  (1 1). 

The  relation  (29)  determines  a  ray  with  a  maximum  or  a  minimum  at  the  height  z^ 
and  entering  the  eye  at  the  height  zq.  If  a  ray  has  a  maximum  or  a  minimum  the  horizontal 
Cartesian  coordinate  x  must  be  a  multi-valued  function  of  the  vertical  coordinate  z,  and,  as 
discussed  in  Section  III  for  the  case  of  a  spherically  symmetric  medium,  continuing  x 
beyond  the  value  given  by  (29)  after  z  reaches  the  value  z^  requires  changing  the 
integration  interval  and  adding  a  constant.  For  programming  the  calculation  in  Mathcad,  a 
somewhat  more  convenient  approach  is  to  solve  a  differential  equation  for  the  ray, 
obtaining  z  as  a  (single  valued)  function  of  x. 


Differentiating  both  sides  of  (29)  provides  a  suitable  differential  equation 

dx  ■  n(zn,) 


(A-2) 


A-1 


The  right  side  of  (A-2)  must  be  positive  if  z  is  increasing  and  negative  if  z  is  decreasing. 
Therefore,  if  the  critical  point  on  the  ray  nearest  the  initial  value  is  a  maximum,^  the  sign 
must  be  positive,  but  if  it  is  a  minimum  the  sign  must  be  negative.  After  integrating  (A-2) 
from  the  eye  position  at  the  point  (0,zo)  to  the  nearest  critical  point,  where  the  right  side 
becomes  infinite,  the  sign  must  change.  If  the  integration  is  continued  the  sign  must 
change  after  each  critical  point  reached  by  the  process. 

Two  versions  of  the  Mathcad  ray  calculation  spread  sheet,  based  on  solving  the 
differential  equation  (A-2),  using  (A-1)  for  the  temperature  profile,  follow.  One,  called 
raycalc3.mcd,  assumes  that  the  critical  point  on  the  ray  nearest  the  eye  level  zq  is  a 
maximum;  the  other,  which  assumes  that  the  nearest  critical  point  is  a  minimum,  is 
raycalc4.mcd.  Both  programs  calculate  Cartesian  coordinates  of  points  along  a  ray  and 
store  the  data  in  a  file.  Each  also  plots  the  corresponding  ray  using  the  data  created  in  this 
manner,  and  a  later  part  of  each  program  can  plot  nine  rays  in  a  single  figure  after  the 
necessary  data  files  have  been  created  by  running  the  earlier  part. 
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The  critical  point  may  occur  at  Zm  or  at  another  value  of  z  at  which  the  index  of  refraction  n(z)  equals 
n(znO. 


A-2 


RAYCALC3.MCD 


Raycalc3.mcd 

CAUTION!  Click  on  Automatic  in  the  Math  menu  to  toggle  to  the  manual  mode  before  proceeding. 

This  program  generates  a  temperature  profile  and  the  corresponding  index  of  refraction  profile.  The 
function  chosen  for  this  purpose  is  a  cubic  polynomial  that  has  only  1  real  zero  and  depends  on  three 
independent  parameters  and  a  scale  factor  sf.  The  program  also  traces  a  ray,  given  its  maximum 
height,  but  it  assumes  a  ray  writh  at  most  one  maximum. 

TOL  =  10  * 

sf  100  b  =  0  c  =  2  d  =  275  a  =  +  c  hi  --  —  sO  =  sf^  a2  =  - 

3  sf^  sf 

t(x)  =  —  -  b2  +  a2  X  +  d  f=—  g  =  2  b2  h=2f 

sf3  sf3 

Check  the  lowest  temperature,  which  is  intitially  in  deg  K.,  in  deg  C  and  deg  F: 

C  :=  t(0)  -  273  C  =  2  9~  +  32.2  =  35.8  C  :=  t(450)  -  273 

9  —  +  32.2  =  167.425 
5 

Get  first  and  second  derivatives  t'(x)  and  t"(x)  of  the  temperature  profile. 
t*(x)  :=  fx^  -  g  X  +  a2  =  h  X  -  g 


t"(z)degKy 


RaycalcS.mcd 


# 


# 


Get  the  index  of  refraction  profile.  First  define  some  universal  constants:  Nq  =  77.5256, 

r  :=  .03417  ,  the  standard  atmospheric  pressure  at  sea  level  Pq  =  1013.25  .which,  along 

with  the  temperature  profile  t(z),  determine  the  atomospheric  pressure  p(z)  and  the  modulus  of 
refraction  N(z)  as  a  function  of  the  height  z.  The  index  of  refraction  n(z)  is  defined  in  terms  of  the 
modified  modulus  M(z).  The  parameter  k  is  .048  if  the  unit  of  height  is  feet  or  .157  if  it  is  meters. 


p(z)  =  Pq  e' 


*z 

-r- 

—  dC 

t(Q 

• 

0 

N(z)  :=No 


P(z) 

t(z) 
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K  :=  .048  M(z)  N(z)  +  K  Z  n(z)  =1  +  10  ''•M(z) 

j  =0,1.  45  Vj  =  10  j  ''"j  ■  ^('j)  u  =  70  X  g  =  root(n(u)  -  n(0)  ,u) 
Write  modified  index  profile  to  a  file.  WRITEPRN(nprofO)  =  auginent(v,w) 

The  first  and  second  derivatives  n'(z),  n"(z)  are  then  given  by: 

*/  \  I  ^  ^  *’(*)  XT/  \\ 

n’(z)  :=  10  [K - ^-^=-^  N(z) 


t(z) 


Ts  =  3T 


n"(z)  =  10' 


2  t’(z)%r3  P(z)-  t(z)  t"(z) 
t(z)^ 


N(z) 


Index  of  Refraction  Profile 


Raycalc3.mcd 

Estimate  the  height  where  the  index  is  a  maximum  and  then  calculate  it.  u  =  200 

h  uj  :=  root(n'(u)  ,u)  h  =  140.33556  Avoid  calculating  a  ray  with  a  maximum  at 

this  height. 

Choose  the  viewing  eye  height:  z  q  =80 

u  :=  400  Ray  maximum  should  be  higher  than  root^n(u)  -  n^ZQ^,u^  =  177.39753 

Calculate  the  height  at  which  the  index  has  the  same  value  as  at  the  ground.  Actually  choose  .3  feet 
above  the  ground  for  this  purpose  to  avoid  complications  in  calculating  the  ray  path. 

u  =  500  h  g  =  root(n(u)  -  n(.3)  ,u)h  g  =  229.76117 

Choose  the  height  ofthe  ray  maximum:  z  ^  :=  300  nuj  :=n^Zjjj^  ”  2  “  ® 

u  =  1  ^0  "  -  "m’“))  ^0  =  ® 


Define  some  auxilliary  functions:  0(z)  I  n(z)  -  ^2 


Calculate  the  ray  arclength  from  the  viewing  eye  to  the  ray  maximum. 


TOL  -  .1 


Z  :=  z  ^  -  1 
mn  m 


n 


m 


mn 


O  ( z 


X 


n 


mn 


m  /  \  1 , 

nfz,„„Vn’(z 


mn 


+  n 


m 


m 


mn 


n”(C) 


(n(C)n’(O) 


XQ=XJ^-X2  *0”  5.75122*  10 

r^i  . 


*00  -  2*0"  "m 


Xqq  =  2.15026-10^ 


<i>(C) 


dC  +  n 


^1  -  ^0-1 


m 


"Ul  -n’U 


n 


m 


0 


00 


^0 


n”(Q 


(n(C)-n’(C)r 


0(C) 


1 


6080 


9.45924 


6080 


35.36618 


Define  the  second  function  used  in  the  ray  differential  equations. 

h(x,z)  :=  if(x<XQ,0(z),-0(z))  g(x,z)  =  if(x<x  oo,h(x,z),<D(z) 


f2(x,z)  :=  if(z<0,0,g(x,z)) 


Differential  equations:  Yq  =  z  q 


D(x,y)  :=  f2(*’yo) 


RaycalcS.mcd 


Solution  of  differential  equations:  Z  =  rkfixed(y  ,0,200000, 200,  D) 

s  =  arclength,  x  =  horizontal  distance,  z  =  vertical  distance 


Ray  Path  for  a  Given  Maximum  Height 


Y  =  augment(x,z)  WRITEPRN(ray8)  =  Y  ^  -  «n 

*0,1 

To  continue  with  another  ray  change  the  ray  filename,  go  to  page  4  and  define  a  new  ray  height. 


To  plot  all  rays  together  first  read  in  data  from  all  ray  files. 


U  =  READPRN(rayl) 

V  =  READPRN(ray2) 

W  =  READPRN(ray3) 

Y  =  READPRN(ray4) 

Z  =  READPRN(ray5) 

P  =  READPRN(ray6) 

Q  -  READPRN(ray7) 

R  =  READPRN(ray8) 

S  =  READPRN(ray9) 

Assign  a  variable  to  the  data  for  each  ray.  v  :  =  200 

Yoi=80 

RaycalcS.mcd 
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L_ 


RAYCALC4.MCD 


# 


Raycalc4.mcd 

CAUTION!  Click  on  Automatic  in  the  Math  menu  to  toggle  to  the  manual  mode  before  proceeding. 


This  program  is  specifically  for  rays  in  an  inferior  mirage,  i.e.,  each  with  a  minimum  and  a  change  in 
the  sign  of  its  curvature.  The  program  generates  a  temperature  profile  and  the  corresponding  index  of 
refreaction  profile.  The  function  chosen  for  this  purpose  is  a  cubic  polynomial  that  has  only  one  real 
root  and  depends  on  three  parameters  and  a  scale  factor  sf.  The  program  also  traces  a  ray,  given  its 
minimum  height,  but  it  assumes  a  ray  with  at  most  one  minimum  and  one  maximum. 


TOL  =  10  * 

y  y  b 

sf  =  100  b  =  0  c  =  2  d  =  275  a  =  b  +  c  b2  = - 

3  3 

t(x)  :=  — —  b2  x^  +  a2  x  -I-  d  f  =  —  g  =  2  b2 

sl3  sl3 


sf3  =  sf^  a2  =  -- 
sf 


h  =  2  f 


Check  the  lowest  temperature,  which  is  intitially  in  deg  K.,  in  deg  C  and  deg  F: 


C  ;=  t(0)  -  273  C  -  2 


32.2  =  35.8 


Get  first  and  second  derivatives  f(x)  and  t"(x)  of  the  temperature  profile. 


t’(x)  :=  fx^  -  g  x  +  a2  t”(*)  “  “  8 
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Raycalc4.mcd 


Get  the  index  of  refraction  profile.  First  define  some  universal  constants:  Nq  =  77.5256, 

r  -  .03417  ,  the  standard  atmospheric  pressure  at  sea  level  Pq  =  1013.25  ,  which,  along 


with  the  temperature  profile  t(z),  determine  the  atomospheric  pressure  p(2)  and  the  modulus  of 
refraction  N(z)  as  a  function  of  the  height  z.  The  index  of  refraction  n(z)  is  defined  in  terms  of  the 
modified  modulus  M(z).  The  parameter  k  is  .048  if  the  unit  of  height  is  feet  or  .157  if  it  is  meters. 


p(z)  =  Pq  e 


•z 

—  dc 

-r- 

t(Q 

« 

0 

K  :=  .048  M(z)  :=  N(z)  +  K  Z 
j  =  0,1.  45  V.  =  10  j 


N(2)  =No-^ 
t(z) 

n(z)  :=  1  +  10"^  M(z) 
w.  :=  MCv.) 

J  \  J  / 


Write  modified  index  profile  to  a  file.  WRITEPRN(  nprofO  )  =  augment  (  v ,  w  ) 
The  first  and  second  derivatives  n'(z),  n"(z)  are  then  given  by; 


n'(z) 

n-(z) 


r  ^  P(z) 
t(z) 


/2t»(z)^  +  r3-t’(z)-t(z)r(z)^ 
\  t(z)^  / 


3r 


•N(z) 


Index  of  Refraction  Profile 
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Raycalc4.mcd 


i  =  0,1.  50  Xj  =  10  i  tj  =  t^x.^  WRITEPRN(temprof)  =  auginent(x  ,t) 
Vj  =  WRITEPRN(nprof)  =  auginent(x,v) 


Estimate  the  height  where  the  index  is  a  maximum  and  then  calculate  it. 
u  :=  200  h  :=  root(n'(u)  ,u)h  ^  =  140.3355589 

Viewing  eye  height:  z  g  =  200  u  =  0  root^n(u)  -  n^ZQ^,uj  =  49.2077951 

Choose  a  minimum  height  for  a  ray:  z  25  "  m  ”  ”(*in)  ®  2  “  ” 

Calculate  the  height  where  the  index  is  the  same, 
u  =  500  ^  g  =  root^n(u)  -  n 

^  g  =  215.6876967  ^  1  =  ^  0  ^ 


Define  some  auxilliary  functions: 

TOL  =  .1 


Calculate  the  ray  arclength  from  the  viewing  eye  to  the  ray  minimum. 
*mn  ^  *in  ^  1 


n 


m 


mn 


1  "in'^(*mn) 

dC  X  2  ■“ 


^(C) 


n 


m 


m 


mn 


n”(C) 


(n(C)n'(O) 


•0(0  dC 


Xo=Xi  +  X2  ^0^  9.8901947*  10 


*00  =  2*0+  "m 


1  .. 

-  d(^ - ; - ; - -y-  +  H 


DCC)  nUi-n'Ui 


m 


0 


n”(C) 


(n(C)n’(O) 


.$(C)dC 


xgg  =  2.1728195-10^ 


*000  -  2  *gg  -  *g 


X  QQQ  “  3.3566196*  10 


Raycalc4.mcd 


Define  the  second  function  used  in  the  ray  differential  equations. 
h(x,z)  if(x<XQ,-(D(z),(5(z))  gl(x,z)  =  if(x>Xoo,-4)(z),h(x,z) 

fl(x,z)  if(x>Xooo,^(z),gi(x,z)) 
f(x,z)  =  if(z<0,0,f  j(x,z)) 

Differential  equations:  -  1q  D(x,y)  =  f  ^x,yQ^ 


Solution  of  differential  equations:  Z  =  rkflxed  ( y ,  0 , 200000 , 200 ,  D  ) 


s  =  arclength,  x  =  horizontal  distance,  z  =  vertical  distance 


k  =  0,1  .200 


6080 


<0.0 


Ray  Path  for  a  Given  Maximum 


Y  :=  augment(x,z)  WRITEPRN(ray9)  =  Y 

To  continue  with  another  ray  change  the  ray  filename,  go  to  page  4,  and  define  a  new  ray  height. 


Raycalc4.mcd 


To  plot  all  rays  together  first  read  in  data  from  all  ray  files. 

U  =  READPRN(rayl)  V  =  READPRN(ray2)  W  =  READPRN  ( ray3 ) 

Y  =  READPRN(ray4)  Z  -  READPRN(ray5)  P  =  READPRN(ray6) 

Q  =  READPRN(ray7)  R  =  READPRN(ray8)  S  =  READPRN(ray9) 


Assign  a  variable  to  the  data  for  each  ray.  v  =200 
i  :=  0,1.  V 


Uf 

=  Uu 

: 

|] 

II 

yi  =  Yj  1 

=  \i 

Pi 

'll 

Ql 

■■i  =  «i,i 

II 

jfi 

o 

II 
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